/// Script to produce relative risk graph and estimates using nlcom command and coefplot package
/// Fabien Cottier, March 2019; updated May 2020
/// SPEI DATA / non-population weight SPEI


///////// Graph immediate SPEI effect only(SPEI Y0)

* write programe for relative risk estimation // margins cannot produce rr estimates for cont. variables
program define rr_est	
	local models ""
	local x=-.5
	
	* Loop to create multiple nlcom statement
	forvalues j=1(1)11 {

		if "`2'"=="linear"{
			local models `"`models' (_b[mean_speiy]*`x')"'
		}
		
		else if "`2'"=="quadratic"{
			local models `"`models' (_b[mean_speiy]*`x'+_b[mean_speiy#mean_speiy]*(`x')^2)"'
		}
		
		local x=`x'+1/10
	}

	* estimates relative risk using nlcom
	est resto `1' 
	local df_est=e(df_r)
	nlcom `models', post df(`df_est') 
	estimates store est_rr
	
	* produce plot using coefplot
	coefplot (est_rr), ///
		eform yscale(range(.75 1.75)) ylabel(.75 "-25%" 1 "0" 1.25 "+25%" 1.5 "+ 50%" 1.75 "+75%") xlabel(1 "-0.5" 3.5 "-0.25" 6 "0" 8.5 "0.25" 11 "0.5") ///
		ytitle(Relative change in irr. migration (%)) xtitle(Standardized Precipitation Evapotranspiration Index) vertical recast(line) ciopts(recast(rarea)) yline(1) 
end




///////// Split sample graph effects SPEI (SPEI Y0)


* write programe for relative risk estimation // margins cannot produce rr estimates for cont. variables
program define rr_est_Sp	// MuL_Sp for multiple lag with split sample
	local models_L0 ""  
	local x=-.5
	* Loop to create multiple nlcom statement
	forvalues j=1(1)11 {
	
		if "`3'"=="linear"{
			local models_L0 `"`models_L0' (_b[mean_speiy]*`x')"'
		}
	
		else if "`3'"=="quadratic"{
			local models_L0 `"`models_L0' (_b[mean_speiy]*`x'+_b[mean_speiy#mean_speiy]*(`x')^2)"'
		}
		
		local x=`x'+1/10
	}
	
	* estimates relative risk using nlcom // With High labor agriculture
	est resto `1'
	local df_est=e(df_r)
	nlcom `models_L0', post df(`df_est')
	estimates store est_rr_L0H

		
	* estimates relative risk using nlcom // With Low labor agriculture
	est resto `2'
	local df_est2=e(df_r)
	nlcom `models_L0', post df(`df_est2')
	estimates store est_rr_L0L
	
	* graph
	coefplot (est_rr_L0H), bylabel(High Agriculture - Immediate effect) || ///
		(est_rr_L0L), bylabel(Low Agriculture - Immediate effect) ||, ///
		eform yscale(range(.75 1.75)) ylabel(.75 "-25%" 1 "0" 1.25 "+25%" 1.5 "+ 50%" 1.75 "+75%") xlabel(1 "-0.5" 3.5 "-0.25" 6 "0" 8.5 "0.25" 11 "0.5") ///
		ytitle(Relative change in irr. migration (%)) xtitle(Standardized Precipitation Evapotranspiration Index) vertical recast(line) ///
		ciopts(recast(rarea)) yline(1)  byopts(row(3))
		
end


